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Abstract 

The spreading of one- and two-component polymer nanodroplets is studied using molecular 
dynamics simulation in a cylindrical geometry. The droplets consist of polymer chains of length 
10, 40, and 100 monomers per chain described by the bead-spring model spreading on a flat 
surface with a surface-coupled Langevin thermostat. Each droplet contains ~ 350000 monomers. 
The dynamics of the individual components of each droplet are analyzed and compared to the 
dynamics of single component droplets for the spreading rates of the precursor foot and bulk 
droplet, the time evolution of the contact angle, and the velocity distribution inside the droplet. 
We derive spreading models for the cylindrical geometry analogous to the kinetic and hydrodynamic 
models previously developed for the spherical geometry and show that hydrodynamic behavior is 
observed at earlier times for the cylindrical geometry. The contact radius is predicted to scale like 
r(t) ~ t 1 / 5 from the kinetic model and r(t) ~ t 1 ^ 7 for the hydrodynamic model in the cylindrical 
geometry. 

PACS numbers: 68.47.Pe 
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I. INTRODUCTION 



Practical applications of the spreading of a liquid on a solid are prevalent in the lubrica- 
tion, coatings, and printing industries, to name a few. Knowledge of the rates of spreading 
and equilibrium configurations of these systems are crucial for improving their performance. 
Extensive experimental, theoretical, and computational work has been undertaken to better 
understand the interaction between a liquid and solid in contact. Most frequently, the liq- 
uids studied are oligomers or polymers in order to remove the influence of evaporation and 
condensation on the droplet spreading dynamics. 

The total energy dissipation in a spreading droplet can be represented as a sum of three 
different components; one due to the hydrodynamic flow in the bulk of the droplet, one due 
to the viscous dissipation in the precursor foot, and one due to the adsorption and desorption 
of molecules to the solid surface in the vicinity of the contact line Experimental measure- 
ments j^l of microscopic droplets compare well with the hydrodynamic model of droplet 
spreading (J, 0, 0, Q] , indicating that hydrodynamic energy dissipation is an important fea- 
ture of droplet spreading. To date, however, simulations fl H 0, H 0, H Q of 

spherical droplets have been unable to approach the droplet size and time duration required 
for hydrodynamic flow to be relevant. 

Although simulations of spreading droplets typically consider a three-dimensional spread- 
ing hemisphere, there are computational advantages for considering a two-dimensional hemi- 
cylinder |17l I18||. The symmetry along the cylinder axis allows periodic boundary conditions 
to be applied in one direction, and a larger droplet radius r can be simulated with fewer 
atoms since the droplet volume scales as r 2 in the cylindrical geometry instead of r 3 for the 
spherical geometry. With larger droplet sizes, this enables us to simulate more viscous sys- 
tems by including polymer chains of length N = 100 for the first time in droplet spreading 
simulations. 

It has been claimed that hydrodynamic dissipation is dominant for small contact angles 
and non-hydro dynamic dissipation is dominant for relatively large contact angles [19] . This is 
reinforced by the fact that for spherical droplets, spreading models have a kinetic dissipation 
term that is linear in the instantaneous contact radius while the hydrodynamic d issip ation 
term has a logarithmic dependence on the instantaneous contact radius We 
show here that for a cylindrical geometry, the hydrodynamic dissipation term is linearly 
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dependent on the contact radius, which suggests that hydrodynamic flow could contribute 
to the dissipation at earlier times for a cylindrical geometry than for the spherical geometry. 
Our simulations show that this is indeed the case. We also show that the r(t) ~ t 1 ^ 10 
scaling of Tanner's spreading law and the r{t) ~ t 1 ^ 7 prediction of molecular- kinetic theory 
for spherical droplets become r{t) ~ t 1 ^ 7 and r(t) ~ t 1 ^ 5 , respectively, in the cylindrical 
geometry. 

Even though many of the liquids used in surface wetting applications are mixtures or sus- 
pensions, most of the research has focused on single component liquids. Some experimental 



|2lL l22l . I23L |2J] and theoretical |2jj |2y] work has been done on binary droplets focusing 
mainly on the equilibrium behavior. Simulations of binary droplets containing from 4 000 
HQ to 25000 monomers have been performed, but larger system sizes are needed to 
adequately model the spreading dynamics. 

In this paper, we present MD simulations of coarse-grained models of one- and two- 
component polymer droplets for chain length N = 10, 40, and 100. These chain lengths 
are chosen since they have a very low vapor pressure and the droplet spreading is not 
influenced by vaporization and condensation. We analyze the dynamics of the components 
of each droplet and compare the spreading rates of two-component droplets to their single- 
component analogues. We derive the equations for the rate of change of contact angle and 
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and hydrodynamic models 



radius for the cylindrical geometry based on kinetic 

4], Q and we use these models to extract physical parameters for each system. 

The paper is organized as follows. Section HH describes the details of the molecular 
dynamics simulations and the application of the Langevin thermostat to the monomers near 
the substrate. It also describes the methods used to analyze the simulation results. Section 
IIHI presents the droplet spreading models for the cylindrical geometry. Section HVl compares 
the spreading behavior of monodispersed droplets in the spherical and cylindrical geometry, 
droplets of different chain lengths, and binary mixtures. The velocity distributions of both 
homogeneous and binary droplets are analyzed in Section and conclusions are presented 
in Section fVTl 
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II. SIMULATION DETAILS 



A. Potentials and Thermostat 

Molecular dynamics (MD) simulations are performed using a coarse-grained model for the 
polymer chains in which the polymer is represented by spherical beads of mass m attached 
by springs. We use a truncated Lennard- Jones (LJ) potential to describe the interaction 
between the monomers. The LJ potential is given by 
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T T 

U LJ (R)={ " ^ " c (1) 

r > r c 

where e and a are the LJ units of energy and length and the cutoff is set to r c = 2.5 o. The 
monomer-monomer interaction e is used as the reference and all monomers have the same 
diameter a. Although in this paper for the binary mixtures we vary only the chain length, 
in future work we will vary the relative interaction strength. For bonded monomers, we 



apply an additional potentia 
elastic (FENE) potential 



where each bond is described by the finite extensible nonlinear 



U F ENE(r) 
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r < Rq 

(2) 
oo r > R 



with k = 30 e and Ro = 1.5 a. 

Droplets composed of polymer chains of length N=10, 40, 100, or 10/40 and 10/100 
mixtures of equal mole fraction of monomers are generated by equilibrating a melt of the 
polymer containing 500 000 monomers at P ~ between two parallel plates in the z direction 
with periodic boundary conditions in the other two directions. The distance between the 
plates L z ~ 90 a. For the cylindrical geometry, the width of the simulation cell in the y 
direction is chosen to be wide enough such that there are no interactions between a chain and 
its periodic image. The larger the width, the better are the statistical averaging of contact 
angle and radius measurements of the droplets. We found that both L y = 40 and 60 a 
give results with reasonable uncertainty in these measurements. For the spherical droplet, 
the dimensions L x = L y . The cylindrical droplets were constructed by removing all chains 
with centers outside of a hemicylinder of radius Ro = 80 a in the xz plane and L y = 40 a, 
which resulted in droplets containing ~ 350 000 monomers. For R = 50 a and L y = 60 a, 



the droplets contained ~ 200 000 monomers. Hemispherical droplets were constructed in a 
similar manner, with initial radii ~ 48 a, resulting in a droplet also of ~ 200 000 monomers. 
The droplets were then placed above a substrate which initially has an interaction strength 
chosen so that the droplet equilibrates with a contact angle near 90°. This is necessary since 
the method of construction of the drop leaves some segments extending into the vapor phase. 
These dangling chain segments quickly coalesce with the droplet after a short equilibration 
run. Hemispheres and hemicylinders were chosen over spheres and cylinders to avoid the 
substantial simulation time required for the isotropic droplet to transition to a cap geometry 
All of the droplets studied here are large enough to avoid the equilibrium contact angle 
dependence on droplet size observed for smaller system sizes 

The substrate is modeled as a flat surface since it was found previously that with the 
proper choice of thermostat, the simulations using a flat surface exhibit the same behavior 
as a realistic atomic substrate. Since simulating a realistic substrate requires several times 
the total number of atoms in the simulation, using the flat surface greatly improves the 
computational efficiency The interactions between the surface and the monomers in the 
droplet at a distance z from the surface are modeled using an integrated LJ potential, 



tHi) 9 -(i) : 



z < z, 



i: 



Ulf\z) = { 6 \- 10 xz/ yz/ J _ " (3) 

z > Zc 

with z c = 2.2a. 

We apply the Langevin thermostat to provide a realistic representation of the transfer of 
energy in the droplet. The Langevin thermostat simulates a heat bath by adding Gaussian 
white noise and friction terms to the equation of motion, 

171^ = -AUi - rriijLTi + Wj(t), (4) 

where m ; is the mass of monomer i, is the friction parameter for the Langevin thermostat, 
— AUi is the force acting on monomer i due to the potentials defined above, and Wj(t) is a 
Gaussian white noise term such that 

(Wi(t) ■ Wjit')) = GksTmaLSijS (t - t') . (5) 

Coupling all of the monomers to the Langevin thermostat would have the unphysical effect 
of screening the hydrodynamic interactions in the droplet and not damping the monomers 



near the surface stronger than those in the bulk. To overcome this, we use a Langevin 
coupling term with a damping rate that decreases exponentially away from the substrate 
35| . We choose the form 



1l{z) =7£exp(<7-z) (6) 

where 7£ is the surface Langevin coupling and z is the distance from the substrate. We 
generally use values of 7^ = 10.0 r _1 and 3.0 r _1 for e w = 2.0 e and 3.0 e, respectively, based 
on earlier work [8] matching the diffusion constant of the precursor foot for flat and atomistic 
substrates. The larger Yl corresponds to an atomistic substrate with larger corrugation and 
hence larger dissipation and slower diffusion near the substrate. 

The equations of motion are integrated using a velo city- Ver let algorithm. We use a time 

1 /9 

step of At = 0.01 r where r = a (—) The simulations are performed at a temperature 



T = e/ks using the lammps code [26j. Most of the simulations were run on 64 to 100 
processors of Sandia's ICC Intel Xeon cluster. One million steps for a wetting drop of 
350 000 monomers takes 24 to 86 hours on 64 processors depending on the radius of the 
droplet. 

B. Analysis details 

For all of the simulations presented here, we extract the instantaneous contact radius 
r(t) and contact angle 9(t) every 400 r. The contact radius is calculated by defining a one- 
dimensional radial distribution function, g{r) = p(r)/p, based on every monomer within 
1.5 a of the surface. The local density at a distance r from the center of mass of the droplet 
is 

P(r) = ^ (7) 
HK J 2LAr v ; 

where N(r) is the number of monomers at a distance between r and r + Ar from the center 
of mass, L is the width of the simulation cell in the periodic direction, and p is the integral 
of p(r) over the entire surface. The contact radius is defined as the distance r at which 
g(r) = 0.98. This approach provides a robust measure of the radius at any point during 
the spreading simulation. The same calculation is used to obtain the droplet radius for ten 
slices of the droplet at incremental heights every 1.5 a from the surface. A line is fit to the 
resulting points and the instantaneous contact angle is determined from the slope of the 




FIG. 1: Molecular weight dependence of the liquid/vapor surface tension for T = e/kg. The solid 
line is a fit to the experimentally observed N~ 2j/3 dependence. 



line. For simulations that exhibit a precursor foot, the monomers within 4.5 a of the surface 
are ignored in the contact angle calculation. 

Fitting the spreading models to the contact angle and contact radius data requires knowl- 
edge of both the surface tension and viscosity. Both involve separate simulations of bulk 
polymers. The surface tension is obtained by simulating the polymer melt in a slab ge- 
ometry so that there are two surfaces perpendicular to the z direction. For the iV > 40 
system, the melt contains 200 000 monomers and each surface has a cross-sectional area of 
4900 cr 2 . For shorter chain lengths, the melts contain 100 000 monomers and each surface 
has a cross- sectional area of 2500 a 2 . After the system equilibrates, the surface tension is 
calculated from the parallel and perpendicular components of the pressure tensor via j^] 



7= 2 



[p± (*) - P\\ 0)] dz. 



(8) 



The driving force for the spontaneous spreading of a droplet is the difference in surface 
tension at each interface. Since the liquid/vapor surface tension 7 is dependent on the chain 
length, the spreading rate is as well. Figure ^ shows a plot of 7 for droplets of several 

observed 



chain lengths obtained from MD simulation. The data fits the experimentally 
molecular weight dependence 7 ~ jV~ 2 / 3 very well and provides a means to both extrapolate 
values of 7 for large molecular weights and estimate the change in spreading rate for different 
chain lengths. 

The surface tensions of the binary droplets are also obtained from mixtures of the two 
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FIG. 2: Surface tension of binary blends of N=40 polymers with N=5 polymers (o) and N=10 
polymers (□) as a function of the bulk mole fraction x^q of monomers on N=40 chains. 



components. To determine the composition dependence of the surface tension, we equilibrate 
blends of N = 40 with N = 5 and with N = 10 at three blend compositions as shown in Fig. 
|2]and N = 100 with iV = 10 at a composition xiqo = 0.5. This allows us to compare the cases 
where there is a large (5/40 system) or a moderate difference (10/40 system) in the surface 
tension of the pure components. For the binary systems, the mixtures contain 200 000 
monomers and each surface has a cross-sectional area of 4900 a 2 . Note that the surface 
tension shown in Fig. El is not a simple mean field (i.e. linear) function of the monomer 
fraction as the fully equilibrated surface composition consists of almost fully shorter chains. 
The viscosity is calculated from the equilibrium fluctuations of the off-diagonal compo- 



nents of the stress tensor 



39j obtained from polymer melt simulations at T = e/ks with the 



bulk pressure P ~ without tail corrections. We do not include the tail corrections to the 
pressure in order to match the system of the spreading droplet. These simulations are run 
up to 84, 000 r. The autocorrelation function of each off-diagonal component of the stress 
tensor is calculated using the Numerical Recipes routine CORREL j^J. The autocorrelation 
functions are av erag ed to improve statistical uncertainty. From this, the viscosity can be 



calculated using 



V 



dt (a a(3 (t)a af3 (0)) 



(9) 



where V is the system volume and cr a p(t) is the a/3 component of the stress tensor at time 
t. The results for 7 and 77 are also summarized in Table HI 

The viscosity of each blend is obtained in the same manner as the pure components. The 
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TABLE I: Bulk properties of bead-spring chains obtained from MD simulation and model fit 
parameters for T = e/ks, P — 0. 



N 


e w (e) 


' Li \ ) 


r \ ) 


7 (e/* 2 ) 


77 (ml to) 


Kinetic 
Co (— ) 

^ u \T<7 J 


Hydro 
a (a) 


Comb 
Co (— ) 


ned 
a (a) 




Xtdro 
s^ftyaro 




10 


2.0 


10.0 


0.869 


0.84±0.02 


11.1±0.4 


72.2 


30.3 


124 


190 


.0014 


.0049 


.0012 


10 


3.0 


3.0 


0.869 


0.84±0.02 


11.1±0.4 


37.4 


60.4 


63.2 


147 


.0029 


.0045 


.0012 


40 


2.0 


10.0 


0.886 


0.94±0.02 


41.7±1.4 


141 


58.2 


339 


181 


.0008 


.0066 


.0011 


40 


3.0 


3.0 


0.886 


0.94±0.02 


41.7±1.4 


59.9 


73.6 


160 


140 


.0037 


.011 


.0018 


100 


2.0 


10.0 


0.892 


0.96±0.02 


132±2 


180 


41.2 


155 


51.8 


.0015 


.0009 


.0009 


100 


3.0 


3.0 


0.892 


0.96±0.02 


132±2 


82.4 


70.7 


417 


124 


.0057 


.013 


.0019 


100 


2.0 


3.0 


0.892 


0.96±0.02 


132±2 


105 


65.0 


678 


155 


.0012 


.016 


.0022 


100 


3.0 


10.0 


0.892 


0.96±0.02 


132±2 


167 


61.2 


126 


77.3 


.0057 


.0005 


.0004 



TABLE II: Bulk properties and model fit parameters for binary droplets. 



N 


e w (e) 


71 (r- 1 ) 


P K 3 ) 


7 (e/<? 2 ) 


r? (m/rc) 


Kinetic 
Co (£) 


Hydro 
a(<j) 


Combi 
Co (^) 


ned 
a(o-) 


Xkin 


Xhydro 


Xcomb 


10/40 


2.0 


10.0 


0.8800 


0.885±0.02 


34.8±1.4 


109 


58.9 


259 


178 


.0009 


.0059 


.0006 


10/40 


3.0 


3.0 


0.8800 


0.885±0.02 


34.8±1.4 


45.5 


73.9 


168 


160 


.0027 


.015 


.0024 


10/100 


3.0 


3.0 


0.8830 


0.90±0.02 


67.2±2.4 


52.8 


75.7 


158 


127 


.0061 


.015 


.0023 



surface tension and viscosity for each blend is given in Table |H] The surface tensions of 
the mixtures are closer to that of the shorter chains since they dominate the liquid/vapor 
interface in the equilibrated system. However, the viscosity of the mixture is more strongly 
influenced by the longer chains. 



III. CYLINDRICAL GEOMETRY DROPLET SPREADING MODELS 

The droplet is modeled as a cylindrical cap as shown by the hatched region in Figure 01 
Here, the droplet volume is defined as the cap of height h and width 2r of the cylinder with 
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FIG. 3: Diagram of the cylindrical cap. 



radius R and length L. The cap height can be expressed in terms of the contact angle, 6, as 



I - cos (0) 
sin (0) ' 



(10) 



Using the definition for the area of a circular segment, A = |i? {29 — sin 20), the radius of 
the cap can be written in terms of and differentiated to give 



dr 
~db 



A 



1/2 



sin 3 



(10 



, cos0 ) — . (11) 

-sin 6> cos 6> J \ 9 - sin 9 cos # / dt 

The free energy is determined by integrating the surface tensions of the liquid/vapor, 
solid/vapor and solid/liquid interfaces over the areas of each interface. For the cylindrical 
cap geometry shown in Fig. El this is given by 



r (t) fdh'ix t)\ 

F {r{t)} = 2r{t)L { 1sl - Isv) + 2 7 L J dx 1 + ( ^ > \ 



1/2 



(12) 



where 7, 'jsv, an d Jsl are the liquid/vapor, solid/vapor, and solid/liquid surface tensions, 
respectively. The height of the cylindrical cap in terms of the cap dimensions is given by 



h'(x, t) 



r(t) 
sin 9 



x 2 sin 2 9 
r(t) 2 



-1/2 



— cos 6* 



(13) 



Combining Eqs. and EI] and differentiating gives 



dF{r(t)} 
dr(t) 



2L 7 



9 



9n 



sin 9 sin 6Y 



(14) 
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Using the standard mechanical description of dissipative system dynamics, the dissipation 
function can be represented as 



dT{r(t);f(t)} _ dF{r(t)} 

*m dr(t) { ' 



where T is the dissipation function jjj, |4l| , which we consider to be composed of a kinetic 
component and a hydrodynamic component T^2 W - The kinetic dissipation term, due to 
molecular adsorption near the contact line, follows the kinetic model introduced by Eyring 
and coworkers [3Q] and applied to spreading of a spherical droplet by Blake and Haynes 
Js]],!^. In the kinetic model, the liquid molecules jump between surface sites separated by 
a distance A with a frequency K. For the spreading cylinder, the velocity of the contact line 
to first order is obtained from Eq. El as 

(16) 

where the friction coefficient Co = ^f T - Here, An is the density of sites on the solid 
surface. Combining Eqs. IT^HTHl we find that the dissipation term due to the surface kinetics 
is 

= Cor(t) 2 L/2 (17) 
The hydrodynamic dissipation term for the spreading droplet is obtained by solving the 



equations of motion and continuity For the spherical droplet, Seaver and Berg |33[ found 
that approximating the spherical cap as a cylindrical disk of the same volume gave results 
that differed from the rigorous derivation by Cox jj| only by insignificant numerical factors. 
We apply the same approximation here, treating the hydrodynamics of the cylinder as 
identical to that of a rectangular box as shown in Fig. 0J Following de Ruijter et al. 
we set the velocity of the upper part of the leading edge to the droplet spreading velocity, 
v x [x = r(t), z = h'} = f(t). With this boundary condition, the velocity profile is simply 

v ' (x ' z) = Wjf t] - (18) 

The hydrodynamic dissipation Ti w is defined as 



V 
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FIG. 4: Rectangular representation of the cylindrical cap. 
Combining Eqs El and EH and integrating gives 



TJ2 = Lr]r{tf (r(t) - a) fti (20) 

W i — i 

where the parameter a has the same meaning as in p] , it is a minimum radius cutoff applied 
to avoid the singularity in the velocity at the z axis. Note that the logarithmic dependence 
of T^2 W on r(t) for the case of a spreading sphere P becomes a linear dependence for the 
case of a spreading cylinder. Since the rectangular box and the cylindrical cap have the 
same volume, we can rewrite Eq. |2U]in terms of the cylinder dimensions 

E, n9 sin 2 9 (r(t) — a) , , 

. = 2i ' r,i) W^4- (21) 

or the spherical geometry, the hydrodynamic dissipation term has been derived previously 



3, 



T = 67rr(t)7/0 [0(t)] r(t) 2 In [r(t) /a] (22) 

where [0(i)] is defined as 

rh\f)(t\] [I + cos 6(t)} sin 6(t) 

TO] " [1- COS 0(f)] [2 + cos 0(f)]' (23) 
We construct a combined kinetic and hydrodynamic model in a manner analogous to de 
Ruijter et al. by combining Eqs. fT4TIT7l and I2T1 



m = 1 (— on 

^ Co 2»(r(t)-a)Bin2fl Uj n g s y I 
2 T r(t)(6»-sin0cos0) x 
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Rewriting this in terms of the contact angle 9 using Eq. HP gives 



d6 (9 — sin6> cos6>\ ( sin 3 # \ 1 T (sine smSo, 

dt = V A J ~ 0-sin0cos0 J Co , wP^jigJ - 

2 ~ r(t)(0-sin | fl£Qa^) 



(27) 



This can be compared to the analogous expressions for a spherical droplet |2j, |3J: 

r(t) = 7 (cos ^o- cos ^) 

Co + 6# [0(t)] In (^) 

d0 / 7T \i/3 (2 - 3 cos + cos 3 0) 4/3 7 (cos 0o- cos 0) 
dt~~ \W) (l-cos0) 2 Co + 6# TO] In 

The kinetic model is obtained by setting ?? = in Eqs. I2"""l through 12*71 and the hydrodynamic 
model is obtained by setting Co = 0. For the kinetic model, the asymptotic solutions of Eqs. 
El and ESI give 

e(t)~(2A) l *(?£) , (28) 



r(t)~2(2A) 2/5 (^) (29) 



6Co. 

1/5 

\6Coj 

as compared to 9(t) ~ t~ 3 / 7 and r(t) ~ t^ 1 ^ 7 for the spherical geometry. Similarly, the 
asymptotic solutions for the hydrodynamic model give 

9(t) ~ (2Af 7 ( 2A , (30) 



48r/ 



r(t)^2(2A) 3 / 7 (^) 1/7 (31) 
as compared to 9(t) ~ t~ 3 / 10 and r(t) ~ t -1 / 10 for the spherical geometry. 

IV. RESULTS 

A. Comparison to spherical geometry 

For wetting droplets, the spreading is characterized by the formation of a precursor foot 
of monolayer thickness that advances ahead of the bulk of the droplet. The bulk region of 
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t (units of t) 



FIG. 5: Radius of the precursor foot for the cylindrical (solid line) and spherical (dashed line) 
geometries and of the bulk droplet for the cylindrical (dotted line) and spherical (dash-dotted 
line) geometries. Both cylindrical and spherical droplets consist of 20 000 polymers of length 
N = 10. The substrate interaction strength e w = 2.0 e is in the fully wetting regime for N = 10, 
7 £ = 10.0 r _1 . 

the droplet follows the precursor foot at a slower spreading rate. This is demonstrated in 
Fig. El where the contact radius of the foot and bulk regions are plotted as a function of time 
for both the cylindrical and spherical geometries. These droplets contain 20 000 chains of 
N = 10 with a substrate interaction strength e w = 2.0e, which is in the fully wetting regime 
for iV = 10. For the otherwise identical systems, the radii of both regions of the cylindrical 
droplet increase faster than those for the spherical droplet. This is a consequence of the 
droplet spreading in one dimension in the cylindrical geometry and two dimensions in the 
spherical geometry. The precursor foot grows diffusively, r 2 (t) ~ t, in both cases. We return 
to further discussion of the time dependence of r{t) in Sec. II V CI 

Figure |B1 shows the time dependence of the contact angle for the same system. For the 
contact angle, the cylindrical and spherical geometries show a comparable spreading rate. 
Since the droplet volume scales as r 3 in the spherical geometry and r 2 in the cylindrical 
geometry, we favor the cylindrical geometry in order to simulate effectively larger droplets 
with the same number of monomers. 

Fits to the kinetic, hydrodynamic, and combined models are performed by taking initial 
guess values for the independent parameters and integrating the expression for d9/dt defined 
in Eqs. EH and |23 for the cylindrical and spherical droplets, respectively. As these data are 
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FIG. 6: (a) Fit of the kinetic (solid line), hydro dynamic (dotted line), and combined (dashed line) 
models to the contact angle data for the cylindrical geometry obtained from MD simulation (o). 
(b) Model fits for the equivalent droplet in the spherical geometry. Both cylindrical and spherical 
droplets consist of 20 000 polymers of length N = 10. The substrate interaction strength e w = 2.0 e 
is in the fully wetting regime for N = 10, 7£ = 10.0 r . 

in the completely wetting regime, the equilibrium contact angle is fixed at 8q = 0°. The 
integration uses the fourth-order Runge-Kutta method to generate a set of data, 6 ca i c (t). The 
parameters are varied using the downhill simplex method until the difference between 
the model and simulation data, \O ca i c (t) — 0(t)\ /0(t), is minimized. The error reported for 
each model is calculated as 



where M is the number of data points in each set of data. 

For the data shown in Fig. |3 the hydrodynamic model provides a more accurate, though 
still only approximate, fit to the data in the cylindrical geometry. The hydrodynamic cutoff 
a ~ 38.1 cr for the spherical geometry and 25.3 a for the cylindrical geometry, indicating 
stronger hydrodynamic dissipation in the cylindrical geometry. For comparison, the friction 
coefficient obtained from the kinetic model Co = 56.3 m/ra for the spherical geometry and 
56.4 m/ra for the cylindrical geometry, are in excellent agreement. 



1 ^ 

* 2 = jfYl 



e C aic(t) - e(t)\ 2 

0{t) 



(32) 
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FIG. 7: Equilibrium contact angle as a function of surface interaction strength for polymer droplets 
composed of N = 10 (o), N = 40 (□), and N = 100 (o) monomers per polymer. 

B. Chain Length Dependence 

The equilibrium contact angles for nonwetting droplets are plotted as a function of the 
surface interaction strength in Fig. [7| The variation of the surface tension with chain length, 
shown in Fig. ^ causes a shift in the wetting transition in terms of the surface interaction 
strength. The contact angles for the N = 10 and N = 40 droplets are taken from earlier work 
l| on spherical droplets. The droplets are large enough to eliminate any equilibrium contact 
angle dependence on the droplet size ^. Contact angles for the iV = 100 droplets are from 
simulations containing 355 000 total monomers. Although the chain length dependence is 
weak for small e w , the wetting transition is shifted to higher e w for larger chain lengths due 
to the increase in the liquid vapor surface tension. The transition occurs near e c w ~ 1.75 e 
for iV = 10 droplets and increases to about e c w ~ 2.25 e for N = 100 droplets. Hence, results 
shown for e w = 3.0 e are in the completely wetting regime for all chain lengths while those 
for e w = 2.0 e are in the completely wetting regime only for N = 10 and N = 40. From 
Fig. the equilibrium contact angle for iV = 100, e w = 2.0 e is 9q ~ 26°. The equilibrium 
contact angle extracted from the kinetic model fit for this droplet is #o — 28°. For droplets 
in the completely wetting regime, the models are fit using an equilibrium contact angle fixed 
at 6 = 0°. 

The fits to the simulation data for various chain lengths for e w = 2.0 e and 7£ = 10.0 r _1 
are shown in Fig. |Sk- Figure |Hb shows results for e w = 3.0 e and 7£ = 3.0 r _1 . The 
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FIG. 8: (a) Fit of the kinetic (solid line), hydrodynamic (dotted line), and combined (dashed line) 
models to the contact angle data for the cylindrical geometry obtained from MD simulation for 
iV = 10 (o), N = 40 (□), and N = 100 (o) droplets for e w = 2.0 e, Y L = 10.0 t" 1 . (b) Same as 
above with e w = 3.0 s, 7£ = 3.0 t" 1 . For clarity, the N = 10 and iV = 40 data sets have been 
shifted upward 60° and 30°, respectively. 

fitting parameters and x 2 values for all of these droplets are listed in Table HJ Overall, the 
combined model produces the best fits to the data, primarily due to the fact that it has two 
fitting parameters while the other two models each have one. The kinetic model also fits 
the data quite well in most cases, which suggests that the combined model overspecifies the 
droplet spreading behavior for these cases. This is reinforced by the fact that the parameters 
extracted from the combined model do not correspond well with the physical system. In 
most cases, the hydrodynamic cutoff a obtained from the combined model is larger than the 
droplet radius and the friction coefficient is larger than that obtained from the kinetic model. 
Previously Q], the single chain diffusion constant was obtained from simulations of polymer 
melts and the friction coefficient was extracted using the Rouse model via D = k R T ' /mNC,R. 
As for the spherical droplets, the friction coefficient obtained from the kinetic model was 
consistently larger than ( R for all cases. 

Although the hydrodynamic model performs better for the cylindrical geometry than 
for the spherical geometry and gives values for a that are less than the droplet radius in 
every case, it still provides the poorest fit to the data of the three models. As seen from 
Fig. |H1 the best fit to the hydrodynamic model is for the system with the highest viscosity, 
N = 100. To explore this point in greater detail, we ran two additional systems, e w = 2.0 e 
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with Yl = 3.0 r _1 and e w = 3.0 e with 7^ = 10.0 r -1 , for A" = 100. The comparison of 
the fits of the three models to all of the iV = 100 systems is shown in Fig. El Only for 
the droplets with the strongest surface dissipation, 7^ = 10.0 r^ 1 , does the hydrodynamic 
model fit the data very well. The hydrodynamic model fit is very poor for 7^ = 3.0 r _1 
regardless of the equilibrium contact angle. The strong surface dissipation slows the surface 
adsorption rate allowing the hydrodynamic behavior to develop in the bulk region of the 
droplet. On an atomic substrate, this is equivalent to increasing the surface corrugation. 
The combined model gives fitting parameters that are comparable to both the kinetic and 
hydrodynamic models only for the two cases where hydrodynamics are important. It should 
be mentioned that the parameters obtained from the combined model are more sensitive to 
the input parameters for the N = 100, 7^ = 10.0 r _1 systems. For these two systems, a 10% 
change in the viscosity results in a 60% change in £0 an d a 15% change in a for the combined 
model. For the other systems, a 10% change in the viscosity results in a 7% change in £0 
and a 1% change in a on average. The kinetic model fits the data better for e w = 2.0 e 
than e w = 3.0 e, presumably because the driving force is smaller since the initial droplet 
is closer to its final equilibrium contact angle. However, unlike the combined model, the 
friction coefficients for these two droplets are consistent with those for the other droplets so 
we consider the fit parameters to be accurate. 

C. Binary droplets 

For binary droplets, the behavior is complicated by the interdiffusion of the two compo- 
nents. As the droplet spreads, the component with the smaller surface tension gradually 
diffuses to the droplet surface. This is evident in the profiles of the binary droplet shown 
in Fig. E3 The droplet, composed of an initial equimolar mixture of monomers belonging 
to chains of length AT = 10 and A" = 40, is on a surface with e w = 2.0 e. Here, the precur- 
sor foot pulls ahead of the bulk region as the droplet spontaneously wets the surface. The 
composition of the precursor foot shows a slight enrichment of the A^ = 10 chains, ranging 
from 59% to 63% in the three frames shown. For e w = 3.0 e, no segregation in the precursor 
foot is observed. This can be understood in terms of the relative distance from the wetting 
transition for each of the components. As shown in Fig. A" = 40 is quite close to the 
wetting transition for e w = 2.0 e and has a significantly slower spreading rate than A" = 10. 



18 




10000 



20000 30000 
t (units of x) 



40000 



FIG. 9: Fit of the kinetic (solid line), hydrodynamic (dotted line), and combined (dashed line) 
models to the contact angle data for the cylindrical geometry obtained from MD simulation for 
N = 100 wither = 2.0 e, 7£ = lCOr" 1 (o), e w = 2.0 e, j s L = 3-Ot- 1 (□), e w = 3.0e, 7£ = lO.Or" 1 
(A), and e w = 3.0 e, 7£ = 3.0 r -1 (o). For clarity, the first three data sets have been shifted upward 
by 60°, 40°, and 20°, respectively. 



However, for e w = 3.0 e, both N = 40 and N = 10 are far from the wetting transition and 
both have fast spreading rates. 

Figure ^2 shows the dynamics of the contact radius of the precursor foot and the bulk 
droplet for the N = 10 and 40 droplets as well as the droplet consisting of an equal monomer 
mole fraction mixture of polymers of chain length N = 10 and 40 for e w = 3.0 e and 
7£ = 3.0 t~ 1 . In general, the spreading rate of the blend falls between that of the two pure 
droplets. For either the foot or the bulk, it does not appear that the dynamics of the blend 
are being dominated by either the N = 10 or iV = 40 polymers. Similar results are found for 
e w = 2.0 e and 7£ = 10.0 r _1 . This is in agreement with earlier simulations 29j,|42j of blends 
of chain length N = 8 and N = 16 where no significant chain length effects are observed. 

The spreading rate of a droplet consisting of an equal monomer mole fraction mixture of 
chains of length iV = 10 and 100 is shown in Fig. El for the case where both components 
are in the completely wetting regime, e w = 3.0 e. Here, the contact radius of the bulk region 
follows the behavior of the iV = 100 droplet more closely than the iV = 10 droplet. This 
indicates that the spreading of the bulk region of the droplet is limited by the diffusion 
rate of the larger component. The spreading rate of the foot is nearly equal to the average 
spreading rate of the two pure components and the precursor foot composition is about 50% 
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FIG. 10: (Color Online) Snapshots of the spreading binary droplet containing 352 000 monomers 
of a mixture of N = 10 and N = 40 polymers. The images are taken from time t = (top), 
t = 40 000 r (middle), and t = 80 000 r (bottom). Monomers from N = 10 chains are shown in red 
and monomers from N = 40 chains are shown in blue. e w = 2.0 e, 7£ = 10. Or -1 . 
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FIG. 11: Spreading rate of (a) the precursor foot and (b) bulk droplet radius for a 352 000 monomer 
mixture of N = 10 and N = 40 polymers compared to homogeneous polymer droplets of the same 
size. The curves correspond to N = 10 (dashed line), N = 40 (dotted line) and the mixture (solid 
line). e w = 3.0 e, 7£ = 3.0 t _1 . 
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FIG. 12: Spreading rate of (a) the precursor foot and (b) bulk droplet radius for a mixture 
of N = 10 and N = 100 polymers compared to homogeneous polymer droplets. The curves 
correspond to N = 10 (dashed line), N = 100 (dotted line) and the mixture (solid line). e w = 3.0 e, 
7 £ = 3.0r- 1 . 

short chain monomers, the same as in the bulk. As found for the 10/40 blend, there is no 
measurable segregation for the more energetic surface. The case where one component wets 
the surface and the other does not will be examined in more detail in another paper where 
we also study the effect of varying e w for two components with the same chain length. 

The three spreading models are fit to the contact angle data for a blend of N = 10 
and N = 40 polymers and for a blend of N = 10 and iV = 100 polymers. The fitting 
parameters and associated errors are given in Table |H] In every case, (o from the kinetic 
model is between that of the corresponding pure component systems. The hydrodynamic 
model gives a value for a that agrees very well with the pure component system of the larger 
chain length. The 10/100 mixture is fit much better by the combined model than by either 
the kinetic or hydrodynamic models, possibly due to a mixture of slow and fast dynamics 
from the two chain lengths. 

V. VELOCITY DISTRIBUTION 

To analyze the droplet spreading dynamics in greater detail, we consider the velocity 
distribution inside the droplet and along the precursor foot. The velocity at a given position 
in the droplet is obtained by generating a histogram of the instantaneous velocities of the 
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FIG. 13: Instantaneous velocity distribution of a homogeneous droplet after 6000 r. N = 40, 
e w = 3.0 e, 7£ = 3.0 r _1 . 

monomers. To eliminate the random fluctuations of the atomic velocities, we average 500 
such histograms over a period of 50 r. After generating the averaged histogram, bins con- 
taining less than 50 monomers are manually removed. Otherwise, these nearly empty bins 
would create the illusion of more flow at the surface than is actually present. 

The instantaneous velocity of a homogeneous droplet composed of iV = 40 polymers is 
shown in Fig. ED The region near the edge of the droplet has the highest velocity while 
the monomers at the center of the droplet and near the substrate but not near the edge 
are almost stationary. This is in sharp contrast to previously published velocity fields of 
spreading droplets |3, llOj which show most of the monomers in the droplet moving at the 
same speed. In these previous cases, the droplets started as equilibrated spheres placed just 
above the substrate, having a contact angle of 180°. Thus, these droplets gain a significant 
amount of momentum as they move towards the substrate and evolve into the cap geometry 
used as the starting point in our simulations. By starting with the droplet equilibrated at a 
contact angle of less than 90°, we are able to focus on the dynamics induced by the surface 
tension driving force and not by the momentum gained by the droplet coming into contact 
with the surface. 

The velocity distribution of each component of an equimolar blend of iV = 10 and iV = 
100 polymers is shown in Fig. El Instead of averaging the instantaneous velocity over a 
short period, the velocity is calculated from the difference in monomer positions at 20 000 r 
and 40 000 r. This more clearly demonstrates the slight differences in spreading behavior for 



22 




x (units of a) 



FIG. 14: (a) Velocity distribution of N = 10 polymers in a 10/100 binary droplet after 20 000 r. 
(b) Velocity distribution of N = 100 polymers in the same droplet. e w = 3.0 e, 7£ = 3.0 t -1 . 

the two chain lengths. Figure El shows that the shorter chains move more rapidly than the 
longer chains, but generally exhibit the same behavior. The shorter chains that are buried 
in the droplet near the substrate show a strong tendency to move away from the substrate 
whereas the longer chains show no such tendency even though the shorter chains are farther 
from the nonwetting transition than the longer chains. This is possibly due to the ability of 
the shorter chains to diffuse a detectable amount during the 20 000 r time period while the 
longer chains are considerably slower. Figure O also indicates that the radial component 
of the velocity increases as the distance from the center of the droplet increases. A more 
detailed analysis of the source material of the precursor foot shows that the radial distance 
traveled by a polymer depends mostly on the initial radial position of the polymer and not 
on how near it is to the droplet surface. Although the highest velocities are found at the 
droplet surface, this does not have much influence on the composition of the precursor foot. 



VI. CONCLUSIONS 



Using molecular dynamics simulation, we study the spreading dynamics of binary polymer 
droplets on a flat substrate. We apply a cylindrical geometry and demonstrate that the 
same qualitative spreading behavior is observed as in a spherical geometry. We derive 
spreading models for the cylindrical geometry based on hydrodynamic and kinetic dissipation 
mechanisms and show that hydrodynamic effects become relevant at shorter time scales in 
the cylindrical geometry. We also show that the r(t) ~ t 1 / 10 scaling from the hydrodynamic 
model and the r(t) ~ t 1 ^ 7 scaling from the kinetic model for spherical droplets become 
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r(t) ~ t 1 / 1 and r(t) ~ t 1//5 , respectively, in the cylindrical geometry. 

Fitting the models to the spreading data of homogeneous droplets shows that the best fit 
is obtained using a combination of the kinetic and hydrodynamic dissipation mechanisms, 
although the parameters extracted from the fit do not agree well with the physical parameters 
of the system. The hydrodynamic model fit the data well only for the slowest spreading and 
highest viscosity (N = 100) case studied, indicating that the kinetic dissipation mechanism 
may be dominating any hydrodynamic effects for the other cases. By increasing the strength 
of the surface dissipation, we are able to slow the droplet spreading rate enough for the 
N = 100 for hydrodynamic dissipation to be significant. 

Compared to homogeneous droplets, the spreading of binary droplets is characterized 
by the difference in surface tension, viscosity, and interaction strength between the two 
components and with the substrate. In the binary droplet at equilibrium, the fraction of 
shorter chains, which have a lower surface tension, is higher at the droplet surface. However, 
the interdiffusion rate is much slower than the spreading rate for the droplets presented here. 
As a result, no enrichment of the lower surface tension component is observed either at the 
droplet surface or in the precursor foot when both components wet the surface. In the case 
that the difference in viscosity between the two components is large, the spreading rate of 
the bulk region is limited by the spreading rate of the more viscous component. Otherwise, 
the spreading rate is roughly equal to the average rate of the two components. The single 
binary system for which the combined model performed noticeably better than either the 
kinetic or hydrodynamic models was the mixture of chain lengths N = 10 and iV = 100, 
possibly due to the combination of fast dynamics from N = 10 chains and slow dynamics 
from N = 100 chains. 

By starting with droplets that have a contact angle 6 = 90° and not with spherical 
droplets above the substrate, we are able to focus on the dynamics induced by the driving 
forces of droplet spreading and not by the momentum gained by coming into contact with 
the substrate. The instantaneous velocity distribution of the spreading droplet shows that 
spreading occurs by the motion of the droplet surface while the interior of the droplet is 
almost stationary. For a droplet composed of an equimolar mixture of short-chain and long- 
chain polymers, we find that the shorter chains move more rapidly than the longer ones near 
the surface of the droplet. From the droplet surface, they move downward to the precursor 
foot and then outward along the substrate. 
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Future work will include studying the effects of the interaction strength between the two 
components and the effects of patterned surfaces. 
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